%% Figure 4: Impact of central clearing on expected default losses with varying loss sharing rules
clear all; clc

%% Create core periphery network and define general help functions
N = 3*5; %number of entities
Create_core_periphery_network
K = 6; %number of derivative classes

pi = 0.05; %default probability
sigma = 1; %idiosyncratic volatility
sigma_M = 0.8; %systematic volatility
beta = 0.2; %systematic risk exposure
alpha_uc = 0.98; %margin confidence level uncleared
alpha_CCP = 0.99; %margin confidence level cleared
xi = @(alpha) (1-alpha).*norminv(1-alpha) + normpdf(norminv(alpha));
f = @(K,beta) sqrt(beta.^2.*sigma_M.^2.*K.^2 + sigma.^2.*K);

%% portfolio characteristics
    eta_K = abs(sum(V,2))./sum(abs(V),2); %net-to-gross ratio
    N_i = sum(abs(V)>0,2); %gross 
    G_K = N_i.*sqrt(beta^2*sigma_M^2+sigma^2); % gross risk
    sigma_bar_K =  N_i .* eta_K .* f(1,beta); %net risk

%% loss sharing rule
    w = @(delta) delta*G_K + (1-delta)*sigma_bar_K; %weights

%% simulate defaults and compute impact of central clearing using Monte Carlo
    nsim = 50000; %number of simulations
    D = rand(1,N,nsim)<pi; %simulated defaults
    
  
    deltas = 0:0.01:1;
    core = NaN*ones(length(deltas),1);
    periphery = core;
    for t=1:length(deltas)
        W = w(deltas(t));
    
        % core 
        i = N/3+1;
        core(t) = f(K-1,beta)/f(K,beta) + W(i)/N_i(i)*xi(alpha_CCP)/xi(alpha_uc)*f(1,beta)/f(K,beta)/pi*mean( sum(D(1,(1:N)~=i,:).*repmat(N_i((1:N)~=i)'.*eta_K((1:N)~=i)',1,1,nsim),2)./(W(i) + sum((1-D(1,(1:N)~=i,:)).*repmat(W((1:N)~=i)',1,1,nsim),2)),3)-1;
        
        % periphery
        i = 1;
        periphery(t) = f(K-1,beta)/f(K,beta) + W(i)/N_i(i)*xi(alpha_CCP)/xi(alpha_uc)*f(1,beta)/f(K,beta)/pi*mean( sum(D(1,(1:N)~=i,:).*repmat(N_i((1:N)~=i)'.*eta_K((1:N)~=i)',1,1,nsim),2)./(W(i) + sum((1-D(1,(1:N)~=i,:)).*repmat(W((1:N)~=i)',1,1,nsim),2)),3)-1;
    end

%% figure
    FontSize = 19;
    linetypes = {'-','-',':','-.','--','--',':'};
    linemarkers = {'none','o','none','.','s','d','+'};

    h1 = figure();
    hold on
    p1=plot(deltas,periphery,'--red','LineWidth',1.2);
    p2=plot(deltas,core,'-black','LineWidth',1.2);
    grid on
    xlabel('Weight of Gross Risk ($\delta$)','FontSize',10,'Interpreter','latex')
    ylabel('$$\Delta$$Expected default loss ($$\Delta DL$$)','FontSize',10,'Interpreter','latex')

    legende=legend([p1,p2],'Periphery','Core','Location','NorthEast','interpreter','latex');
    set(legende,'FontSize',FontSize,'Interpreter','latex','color','none');
    set(gca,'FontSize',round(FontSize*4/5));

    set(h1,'PaperPositionMode','auto');
    set(h1,'Position',[100 100 600 400]);

    ylabels=get(gca,'YTickLabel'); % get the y axis labels
    for i=1:size(ylabels,1)
        ylabels_modif{i}=sprintf('%.3g%%', str2double(ylabels{i})*100);
    end
    set(gca,'YTickLabel',ylabels_modif)

    print(h1, '-depsc2','-loose','pics/Figure_4.eps')